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Abstract 

Iterative  methods  for  the  linear  systems  of  algebraic  equations  arising 
from  elliptic  finite  element  problems  are  considered.  Methods  previously 
known  to  work  well  for  positive  definite,  symmetric  problems  are  extended 
to  certain  nonsymmetric  problems,  which  also  can  have  some  eigenvalues 
in  the  left  half  plane. 

We  first  consider  an  additive  Schwarz  method  applied  to  linear,  sec- 
ond order,  symmetric  or  nonsymmetric,  indefinite  elliptic  boundary  value 
problems  in  two  and  three  dimensions.  An  alternative  linear  system,  which 
has  the  same  solution  as  the  original  problem,  is  derived  and  this  system 
is  then  solved  by  using  GMRES,  an  iterative  method  of  conjugate  gradient 
type.  In  each  iteration  step,  a  coarse  mesh  finite  element  problem  and  a 
number  of  local  problems  are  solved  on  small,  overlapping  subregions  into 
which  the  original  region  is  subdivided.  We  show  that  the  rate  of  conver- 
gence is  independent  of  the  number  of  degrees  of  freedom  and  the  number 
of  local  problems  if  the  coarse  mesh  is  fine  enough.  The  performance  of 
the  method  is  illustrated  by  results  of  several  numerical  experiments. 

We  also  consider  two  other  iterative  method  for  solving  the  same  class 
of  elliptic  problems  in  two  dimensions.  Using  an  observation  of  Dryja  and 
Widlund,  we  show  that  the  rate  of  convergence  of  certain  iterative  sub- 
structuring  methods  deteriorates  only  quite  slowly  when  the  local  problems 
increase  in  size.  A  similar  result  is  established  for  Yserentant's  hierarchical 
basis  method. 

Key  words  Schwarz's  alternating  method,  domain  decomposition,  nonsymmetric 
and  indefinite,  elliptic  equations,  finite  elements 
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1      Introduction 

Domain  decomposition  techniques  are  powerful  iterative  methods  for  solving  lin- 
ear systems  of  equations  that  arise  from  finite  element  problems.  In  each  iteration 
step,  a  coarse  mesh  finite  element  problem  and  a  number  of  smaller  hnear  sys- 
tems, which  correspond  to  the  restriction  of  the  original  problem  to  subregions, 
are  solved  instead  of  the  large  original  system.  These  algorithms  can  be  regarded 
as  divide  and  conquer  methods.  The  number  of  subproblems  can  be  large  and 
these  methods  are  therefore  promising  for  parallel  computation.  The  central 
mathematical  question  is  to  obtain  estimates  on  the  rate  of  convergence  of  the 
iteration  by  deriving  bounds  on  the  spectrum  of  the  iteration  operator.  We  are 
able  to  establish  quite  satisfactory  bounds  if  the  coarse  mesh  is  fine  enough. 

We  work  with  two  triangulations  of  the  region  into  subregions,  also  called 
substructures,  which  define  a  coarse,  global  model,  and  into  elements  of  a  finite 
element  model.  As  in  the  positive  definite  case  considered  previously,  see  Cai 
[3], [4],  Dryja  [6]  and  Dryja  and  Widlund  [7], [8],  the  coarse  problem  provides 
interchange  of  information  between  the  different  parts  of  the  region;  it  is  known 
that  without  such  a  coarse  subproblem  the  rate  of  convergence  is  considerably 
much  slower.  This  part  of  the  approximate  solver  plays  an  additional  role  in 
the  indefinite  case.  We  can  interpret  the  main  resiilts  of  this  paper  by  saying 
that  if  the  eigenfunctions  corresponding  to  the  eigenvalues  in  the  left  half  plane 
are  approximated  well  enough  on  the  coarse  mesh,  then  the  spectrum  of  the 
preconditioned  linear  system  of  equations  lies  in  a  fixed  bounded  subset  of  the 
right  half  plane.  This  is  important  for  the  rate  of  convergence  of  the  iterative 
method.  The  least  favorable  situation  for  iterative  methods  of  conjugate  gradient 
type  is  the  case  where  the  origin  of  the  complex  plane  is  surrounded  by  eigen^•alues 
of  the  iteration  operator.  Here  we  are  able  to  avoid  such  a  situation. 

The  additive  Schwarz  algorithms,  introduced  in  [7],  cf.  also  [6], [8], [9], [18], 
provide  a  means  of  constructing  preconditioners  for  many  problems  in  terms  of 
a  partition  of  a  given  finite  element  space  into  a  sum  of  subspaces.  The  use  of 
such  a  preconditioner  involves  solving,  exactly  or  approximately,  the  restriction 
of  the  original  problem  to  the  different  subspaces.  The  residual,  which  plays  a 
central  role  in  the  iteration,  is  computed  as  a  sum  of  terms  from  the  different 
subspaces.  These  terms  can  be  computed  in  parallel.  We  note  that  it  has  been 
shown  in  Dryja  and  Widlund  [9]  that  many  domain  decomposition  methods  can 
be  viewed  as  additive  Schwarz  methods.  For  recent  work  on  the  case  of  more 
than  two  levels  of  triangulation,  see  Dryja  and  Widlund  [10]  and  Xu  [23]. 

In  the  symmetric,  positive  definite  case,  the  iterative  method  most  commoiUy 
used  to  solve  the  transformed  (preconditioned)  equations  is  the  conjugate  gra- 
dient method.  For  the  cases  considered  here  symmetry  is  always  lost.  In  our 
experiments,  we  have  used  a  generalized  conjugate  residual  method  GMRES. 
Since  the  spectrum  of  the  operator  is  confined  to  the  right  half  plane,  Manteuf- 


fel's  Chebyshev  algorithm  would  also  be  successful;  cf.  [17].  Since  we  can  show 
that  the  symmetric  part  of  the  operator  is  iiniformly  positive  definite,  with  re- 
spect to  a  suitable  inner  product,  and  that  the  spectrum  is  uniformly  bounded, 
we  can  guarantee  a  rate  of  convergence,  which  is  independent  of  the  mesh  size 
and  the  number  of  subregions. 

Other  methods  for  indefinite,  eUiptic  problems  are  discxissed  in  [2], [14], [16]. 

The  paper  is  organized  as  follows.  In  Section  2,  we  introduce  the  indefinite, 
eUiptic  boundary  value  problem,  the  two  triangtilations  of  the  domain  and  a 
Galerkin  finite  element  method.  We  briefly  review  the  GMRES  method  in  Section 
3.  In  Section  4,  we  present  two  variants  of  the  additive  Schwarz  method  and  a 
detciiled  analysis  of  their  rates  of  convergence.  Our  ajialysis  is  based  on  previous 
work  on  the  positive  definite  case  and  a  result  due  to  Schatz  [22].  Schatz's  work, 
in  turn,  is  based  on  Garding's  inequaUty  and  the  Aubin-Nitsche  trick;  see  Ciarlet 
[5]  or  Nitsche  [20].  In  Section  5,  we  discuss  some  numerical  results.  Finally,  in 
Section  6,  we  show  that,  for  problems  in  the  plane,  our  result  can  be  extended  to 
iterative  substructuring  and  hierarchical  basis  algorithms  discussed  in  Dryja  and 
Widlund  [8], [9]  and  Yserentant  [24],  respectively. 

2      The  Elliptic  Problems 

Let  n  be  an  open,  bounded  polygonal  region  in  i?'',  tf  =  2  or  3,  with  boundary 
dCl.  Consider  the  homogeneous  Dirichlet  boundary  value  problem: 

{Lu    =    f    in     n, 
u    =    0     on    dn.  ^^' 

The  elliptic  operator  L  has  the  form 

i-C-)  =  -  E  |-(Mx)^)  +  2I:Mx)^  +  c(x)u(x). 

All  the  coefficients  are,  by  assumption,  sufficiently  smooth  and  the  matrix 
{a,_,(i)}  is  symmetric  and  uniformly  positive  definite  for  Vi  6  17.  The  right  hand 
side  /  e  L^{n).  We  also  assume  that  the  equation  haiS  a  unique  solution  in  Hl{Q,). 

Let  (•,  •)  denote  the  usual  L^  inner  product  and  ||  •  ||  or  ||  •  ||l2  the  corresponding 
norm.  The  weak  form  of  equation  (1)  is:  Find  u  e  HKil)  such  that 

B{u,v)  =  {f,v),       \/veH'o{n).  (2) 

The  bilinear  form  B{u,  v)  is  defined  by 

B{u,v)  =2^        Q.Jo      o    dx  +  }22      bi--vdx+  /    cuvdx 
i~^i  J  n      oxj  dxi  ^    J(i     dx,  J  n 


„,       ,         ^    /■         du  dv  ^        ^  f  ^  du         d{b.u)     ,      ,    f   . 


Here,  c{x)  =  c(x)  -  Ef=i  db.{x)/dxi. 

We  also  use  two  other  bilinear  forms 

JU    r         du  dv 


Mu,v)  =    j:J^a,,£.£-dx 


and 


S{ 


.        ^  r  ,  du        dibiu)    . 
u,v)  =   2^  /  o,-—v+  — vdx 

~[  Jil       ox,  OXi 


which  correspond  to  the  second  order  terms  and  the  skew-symmetric  part  of  L, 
respectively.  The  biUnear  form  A  defines  a  norm,  which  we  denote  by  ||  •  ||  a- 
Under  the  assumptions  on  the  coefficients  Qjj,  this  norm  is  equivalent  to  the  Hq 
norm.  It  is  also  easy  to  verify  that 

s(u,r)  =  -s{v,u),  yu,veHlin). 

Throughout  this  paper,  c  and  C,  with  or  without  subscripts,  denote  generic, 
strictly  positive  constants.  They  are  independent  of  the  mesh  parameters,  which 
will  be  introduced  later  in  this  section. 

Using  elementary,  standard  tools,  it  is  easy  to  establish  the  following  inequal- 
ities: 

(i)  I  B{u,v)  \<  C||ulU||t;|U,    Vu,t;  G  HH^). 

(ii)  Garding's  inequality:  There  exists  a  constant  C,  such  that 

ll"ir^  -  C|l"irL.(n)  <  B{u,u),   Vu  G  HUil)  . 

(iii)  There  exists  a  constant  C,  such  that 

I  S{u,v)  \<  C\\u\\A\\v\\i^^n),   Vu,r  6  Hlin), 
I  S(u,t;)  |<  C||t;|U||ullz^(n),    Vu,t;  G  HHil). 

We  note  that  the  bounds  for  B{-,  •)  and  S(-,  •)  are  different,  since  each  of  the 
terms  in  5(-,  •)  contains  a  factor,  which  is  of  zero  order.  This  enables  us  to  control 
the  skew-symmetric  term  and  makes  our  analysis  possible. 


We  also  use  the  following  regularity  resiilt;  cf.  Grisvard  [13]  and  Necas  [19]. 
(iv)  The  solution  w  of 

B{cf>,w)  =  (p,</.),  v<^G//^(n) 

satisfies 

||u;||//>+^(n)  <  C||p||L2(n)  , 

where  7  depends  on  the  interior  angles  of  dO.,  is  independent  of  g  and  is  at  least 

1/2. 

We  approximate  equation  (2)  by  a  Galerkin  conforming  finite  element  method. 
For  simplicity,  we  consider  only  continuous,  piecewise  linear,  triangular  elements 
in  R^  and  tetrahedral  elements  in  R^. 

To  define  the  additive  Schwarz  algorithms,  we  need  two  levels  of  triangulation 
that  have  already  been  introduced  in  [3],  [4], [6], [7], [8], [9].  We  first  partition  Cl 
into  substructures  {^i}-,  which  provide  a  regular  finite  element  triangulation  of 
n.  The  f2,  are  non-overlapping,  <f-dimensional  simpUces.  They  satisfy  all  the 
standard  rules  of  finite  elements;  cf.  Ciarlet  [5].  This  is  the  coarse  mesh  and  it 
defines  a  mesh  parameter  if  =  max{.fifi,  •••,  H;v}- The  triangulation  is  assumed 
to  be  shape  regular,  i.e.  if,,  the  diameter  of  n,  is  bounded  uniformly  in  terms  of 
the  diameter  of  the  largest  inscribed  ball  in  fi,. 

In  a  second  step,  we  divide  each  substructure  fi,  into  smaller  simpUces,  de- 
noted by  {  t/,  j  =  I,---  }.  They  form  a  shape  regular,  fine  mesh  (/i-level)  finite 
element  triangulation  of  fi  with  the  mesh  parameter  h  —  max,,j{/i;-}.  Here  h\'\% 
the  diameter  of  r/. 

We  can  now  define  the  piecewise  linear  finite  element  spaces  over  the  if-level 
and  the  /i-levd  triangulations  of  fl. 

V^  =  {v^  I   continuous  on  fl,  v^  |n,   linear  ,  v^  =  0  on  dD.)  , 

V''  =  {u''  I    continuous  on  fi,  t;'*  [^    linear  ,  u''  =  0  on  5fl}  . 

The  Galerkin  approximation  of  equation  (2)  is  defined  by:  Find  u**  6  V*, 
such  that 

B(u\t;'')   =   i/y),    V^;''GV^  (3) 

If  the  mesh  size  h  is  small  enough,  it  follows  from  a  result  by  Schatz  [22]  that 
this  problem  has  a  unique  solution.  By  using  a  nodal  basis  of  the  finite  element 
space,  equation  (3)  is  transformed  into  a  linear  system  of  algebraic  equations, 
which  is  large,  sparse,  nonsymmetric,  indefinite  and  relatively  ill-conditioned. 


3       A  Brief  Discussion  of  the  GMRES  Method 

Among  the  possible  iterative  methods  to  solve  the  linear  system,  we  have  only 
used  one,  the  GMRES  method;  cf.  Saad  and  Schultz  [21]  and  Eisenstat,  Elman 
and  Schultz  [11].  This  is  a  generalized  minimum  residual  method,  which  in  prac- 
tice has  proven  quite  powerful  for  a  large  class  of  nonsymmetric  problems.  The 
GMRES  method  is  described  in  [21]  and  the  theory  developed  in  L^{9.)  can  be 
found  in  [11].  Both  the  algorithm  and  the  theory  can  easily  be  extended  to  an 
arbitrary  Hilbert  space.  We  briefly  describe  the  algorithm  and  state  a  theorem 
without  proof. 

Let  P  be  a  linear  operator  in  the  finite  dimensional  space  R"  with  an  inner 
product  [•,•],  and  a  corresponding  norm  ||  •  ||,  chosen  to  take  advantage  of  the 
special  properties  of  P.  In  our  appHcations,  P  is  the  preconditioned  stiffness 
matrix  and  the  A-norm  is  used.  P  is  not  symmetric  but  is  positive  definite 
with  respect  to  [•,•].  The  GMRES  method  is  used  to  solve  the  linear  system  of 
equations 

Px   =   6, 

where  6  G  /?"  is  given.  We  begin  from  an  initial  approximation  xq  6  P"  and  the 
initial  residual  tq  =  b  —  Pxq.  In  the  m""  iteration,  a  correction  vector  Zm  is 
computed  in  the  Krylov  subspace 

ICmiro)   =   span{ro,Pro,---,P'"-Vo}, 

which  minimizes  the  norm  of  the  residual.  In  other  words,  Zm  solves 

min     116  -  P(xo  +  ^)||. 

i6ACm('-o) 

The  77i"*  iterate  is  x^   =   Xo  +  Zm- 

The  exact  solution  would  be  reached  in  no  more  than  n  iterations  if  we  use 
exact  arithmetic. 

FoUowing  Eisenstat,  Elman  and  Schultz  [11],  the  rate  of  convergence  of  the 
GMRES  method  can  be  characterized  in  terms  of  the  minimal  eigenvalue  of  the 
symmetric  part  of  the  operator  and  the  norm  of  the  operator.  They  are  defined 

by 

Cp  -   mf  -. T-    and     C^   =  sup    ..    ,.   ■ 

^         x,£o    [x,x]  x^io    Ijxlj 

By  considering  the  decrease  of  the  norm  of  the  residual  in  a  single  step,  the 
following  theorem  can  be  established. 

Theorem  ('E'tsensiai,  Elman  and  Schultz).  If  Cp  >  0,  then,  the  GMRES 
method  converges  and  after  m  steps,  the  norm  of  the  residual  is  bounded  by 

lir.il    <    (1   -    c&)        Ikoll- 


4      Algorithms  on  Overlapping  Subregions 

In  this  section,  we  introduce  two  variants  of  an  additive  Schwarz  algorithm  and 
provide  bounds  on  their  convergence  rates;  see  Theorem  1.  The  analysis  is  valid 
for  both  two  and  three  dimensions. 

We  first  form  a  basic  decomposition  of  the  domain  Q,  into  overlapping  sub- 
regions  and  then  introduce  the  projections  in  terms  of  which  our  algorithms  are 
defined. 

We  use  the  ff -level  subdivision  {Jl,}  of  fi.  Each  subregion  ft,  is  extended  to 
a  larger  region  $7-,  i.e.  ft,  C  ft',.  The  overlap  is  generous  in  the  sense  that  there 
exists  a  constant  a  >  0,  such  that 

distance(an  ■  n  Q,  dfl,  n  fi)  >   aH„    Vz. 

We  assume  that  5fi-  does  not  cut  through  any  ^-level  elements.  We  use  the 
same  construction  for  the  subregions  that  intersect  the  boundary  dCl  except  that 
we  cut  off  the  part  that  is  outside  Q. 

We  also  use  the  notation  fig  =   ^• 

We  note  that  the  larger  the  a,  the  fewer  iterations  can  be  expected.  However, 
if  we  increase  the  overlap,  the  size  of  the  subproblems  increases  and  with  it  the 
cost  of  the  subproblems.  It  is  an  important  practical  issue  to  balance  the  total 
number  of  iterations  and  the  cost  of  solving  the  subproblems. 

For  each  fi,,  a  regTilar  finite  element  subdivision  is  inherited  from  the  /^-level 
subdivision  of  d.  The  corresponding  finite  element  space  is  defined  by 

The  elements  of  this  subspace  of  V*  can  be  extended  continuously  by  zero  to  the 
complement  of  n[.  We  also  use  the  subspace 

Vo''  =  V". 

It  is  easy  to  see  that  our  finite  element  function  space  V*  can  be  represented 
as  the  sum  of  the  A''  -|- 1  subspaces, 

y'>  =  v^^  +  y^''  +...+  v^. 

We  can  now  define  the  projection  operators,  which  are  the  main  building 
blocks  of  our  algorithms.  These  operators  map  the  finite  element  space  V*  onto 
the  subspaces  V/*  and  are  defined  in  terms  of  the  bilinear  forms  B{;  ■)  and  A{-,  •). 


Definition:     Fori  =  0,--',N: 

For  any  w^  6  V*,  QiW^  G  V/*  is  the  solution  of  the  finite  element  equation 

B{Qiw\v!:)   =   B{w\v'^),    Vt;feV;\ 

For  any  w'*  €  V",  PiW^  G  V/"  is  the  solution  of  the  finite  element  equation 

We  now  introduce  the  two  operators,  which  define  our  transformed  equations 

and 

g(2)  =  Qo  +  Pi  +...  +  p^. 

Our  mzdn  effort  goes  into  the  study  of  the  spectra  of  these  two  operators.  The 
only  difference  between  Q^^^  and  Q*^^  is  that,  for  i  >  0,  we  replace  the  projection 
Q,,  corresponding  to  fi,,  by  P,.  The  coarse  mesh  projection  is  not  changed. 

The  computation  of  Qiw'*  or  PiW^^  for  t  >  0  and  an  arbitrary  function 
w^  G  V*,  involves  the  solution  of  a  standard  finite  element  linear  system  of 
algebraic  equations  on  the  small  subregion  Q.-.  The  former  gives  rise  to  a  nonsym- 
metric  linear  system  of  equations  and  the  latter  to  a  positive  definite,  symmetric 
problem.  For  i  =  0,  the  problem  is  a  standard  finite  element  equation  on  the  H- 
level,  coarse  space.  One  can  view  P^  as  a  preconditioner  of  Qi  in  the  subspace  V/*; 
cf.  the  discussion  in  Dryja  and  Widlund  [8], [9].  The  cost  of  the  computation  can 
often  be  decreased  by  introducing  such  simplified  problems.  In  particular,  if  it  is 
possible  to  choose  some  of  the  n[  to  be  rectangular  and  the  corresponding  mesh 
to  be  uniform,  a  Fast  Poisson  solver  can  be  used  to  compute  the  contribution 
from  V/*.  It  is  an  easy  exercise  to  modify  our  theory  to  cover  such  a  case. 

We  will  consider  two  additive  Schwaxz  algorithms: 

Algorithm  1:  Obtain  the  solution  of  equation  (S)  by  solving  the  equation 

Q(^)u'*  =  b^'\  (4) 

and 

Algorithm  2:  Obtain  the  solution  of  equation  (S)  by  solving  the  equation 

Q(2)u''  =  b^'\  (5) 


In  order  for  equations  (4)  and  (5)  to  have  a  unique  solution,  the  operators 
Q^^^  and  Q'^)  must  be  invertible.  This  follows  from  Theorem  1.  To  obtain  the 
same  solution  as  equation  (3),  the  right  hand  sides  6*^^  and  6'^^  must  be  chosen 
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correctly.  The  crucial  observation  is  that  these  right  hand  sides  can  be  computed 
without  knowledge  of  the  solution  of  equation  (3).  The  following  formulas  are 
valid: 

6(1)    =   QW^h    =   J2q,u'^ 

•=0 

and 

6(2)   =   qWu'  =   Qou'  +  J2P'^'- 
«=i 
Each  of  these  terms  can  be  computed  by  solving  a  problem  in  a  subspace,  since 
by  equation  (3), 

and 

A{P.u\v'^)  =  B{u\v^)  =  {f,v^),    Vi^fel^A 

The  main  result  of  this  paper  is  Theorem  1.  By  combining  it  with  the  Theorem 
given  in  Section  3,  we  establish  that  the  two  algorithms  converge  at  a  rate  which 
is  independent  of  the  mesh  parameters  h  and  H,  if  the  coarse  mesh  is  fine  enough. 

Theorem  1    There  exist  constants  Ho  >  0,  c{Ho)    >    0  and  C{Ho),  such  that  if 
H  <  Ho,  then,  for  i  =  1,2 

and 

A(Q('V,(3('\'')  <  c(/fo)A(u^u'•). 

The  special  constant  Co  is  introduced  in  Lemma  1. 
Remarks: 

(a)  The  operator  Qo  is  very  important,  since  it  provides  global  transportation 
of  information.  All  the  other  projections  are  local  mappings.  Without  using 
Qo,  information  would  travel  only  from  one  subregion  to  its  neighbors  in  each 
iteration  and  it  would  take  0{1/H)  iterations  for  the  information  to  propagate 
across  the  region.  Without  such  a  global  mechanism,  it  would  also  be  impossible 
to  confine  the  spectrum  to  the  right  half  plane. 

(b)  The  constant  Hq  determines  the  minimal  size  of  the  coarse  mesh  problem 
and  it  depends  on  the  operator  L.  In  general.  Ho  decreases  if  we  increase  the 
coefficients  of  the  skew-symmetric  terms,  it  decreases  with  c,  while  it  increases 
if  we  increase  the  overlap.  Hq  also  depends  on  the  shape  of  the  domain  fi.  If 
the  domain  is  not  convex,  the  estimate  of  Hq,  imphcit  in  our  proof  of  Lemma  5, 
depends  on  the  parameter  7.  We  do  not  have  an  explicit  formula  for  Ho  but  we 
know  from  experience  that  it  can  be  determined  by  numerical  experiments. 

(c)  K  the  operator  L  is  positive  definite,  synmietric,  there  is  no  restriction  on 
the  coarse  mesh  size  H,  i.e.  Ho  =  cx). 

The  proof  of  Theorem  1  is  based  on  the  following  results. 


Lemma  1    There  exists  a  constant  Cq,  which  is  independent  of  h  and  H,  such 
that,  for  all  u''  e  V*,  there  exist  uf  G  V;'*  xvith 

N 
1=0 

and 

»=0 

This  lemma  is  also  central  in  the  theorj'  previously  developed  for  positive 
definite,  symmetric  problems.  For  a  proof  see  Dryja  and  Widlund  [8];  cf.  also 
Lions  [15].  Note  that  this  lemma  is  independent  of  the  skew-symmetric  and  zero 
order  terms  of  the  elliptic  operator.  In  the  symmetric,  positive  definite  case, 
Lemma  1  is  combined  with  an  abstract  argument  to  give  a  lower  bound  for  the 
spectrum  of  the  iteration  operator. 

The  next  lemma  is  a  variation  of  a  result  by  Schatz;  cf.  [22].  In  his  proof, 
Garding's  inequality,  (ii),  and  the  regularity  result,  (iv),  are  used.  The  proof 
of  Lemma  2  follows  directly  from  Schatz's  work  by  replacing  the  approximate 
solution  by  the  coarse  mesh  solution  and  the  exact  solution  of  the  continuous 
problem  by  the  finite  element  solution  in  V^. 

Lemma  2    There  exist  constants  Ho  >  0  and  C{Ho),  such  that  if  H  <  Hq,  then, 

HQou'^Wa  <  C{Ho)\\u^\\a 

and 

IIQou''  -  u'^Wl^  <  C{Ho)H''\\Qou''  -  u'^Wa- 

Lemma  3   The  restriction  of  the  quadratic  form  B{-,  •)  to  the  subspaces  V-^,  i  > 
0,  is  strictly  positive  definite  for  H  sufficiently  small,  i.e. 

cA{u\u'')  <  B(u\u''),Vu''  e  V^  . 

Proof  of  Lemma  3:  We  have  to  prove  that  the  second  order  term  dominate 
the  other  symmetric  term.  This  follows  from  the  fact  that  the  smallest  eigenvalue 
for  the  Dirichlet  problem  for  —A  on  the  region  Q,[  is  on  the  order  of  H~^. 

Lemma  4  Let  v'*  =  J^v^  ,  where  v^  €  V/*.  Then  there  exists  a  constant  C,  such 
that 


•'fA- 


Proof  of  Lemma  4:  The  proof  follows  from  the  observation  that  for  each 
X  6  n,  the  number  of  terms  in  the  sum,  which  differ  from  zero,  is  uniformly 
bounded. 
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Lemma  5  There  exist  constants  Hq  >  0,  c{Ho)  and  C{Hq),  such  that  if  H  <  Hq, 
then, 

N 

c{Ho)Co'A{u\u'')  <  J2MQ>^\Qi^'')  <  C{Ho)A{u\u'') 

1=0 

and 

c{Ho)Co^A{u\ u*")    <    AiQou\ Qou'*)  +  Eili  A{Piu\ P.«'') 

<  C{Hq)A{u'',u''). 

Proof  of  Lemma  5:  An  upper  bound  for  A{Qou'^,  Qou'*)  is  given  in  Lemma 
2.  To  obtain  an  upper  bound  for  the  sum  of  the  other  terms,  we  iise  Lemma  3 
and 

B(QiU^Q.u'')  =  B(u^Q.u'') 

to  show  that 

cf2A{Q.u\Qiu'')  <  B{u\f2Q,u'). 

The  right  hand  side  can  be  estimated  by  using  inequality  (i)  and  Lemma  4.  The 
other  upper  bound  is  established  in  a  similar  way. 

To  prove  the  lower  bounds,  we  begin  by  using  Lemma  2  and  the  triangle 
inequality  to  obtain 

Wu'W]^  <  C{H'-'A{u\u')  +  WQou'Wh). 

Since  the  eigenvalues  of  the  Dirichlet  problem  for  —A  are  bounded  from  below 
and  Lemma  2  holds,  the  last  term  can  be  replaced  by  C||(5ou''|U||"''|U-  By  using 
Garding's  inequality,  (ii),  it  follows  that 

(1  -  CH'-')A{u\u'')  <  B{u\u'')  +  C||(?ou''|U||u''|U  • 

By  the  definition  of  the  operators  Q,  and  Lemma  1,  we  find  that 

B{u\  ti")  =  Yl  B{u\  ul)  =  X:  B{Q,u\  uj-). 

t=0  »=0 

The  boundedness  of  B{-,  •),  (i),  can  now  be  used  to  obtain 

f:B(Q,u\u?)<cf:iic?.«''iuiiufiu, 

«=0  «=0 

which  by  Lemma  1  can  be  estimated  by 

t=0 
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We  finally  obtain 


Aiu^u'^)   <   CCl  f^A{Q,u\Q,u''), 

i=0 


for  sufficiently  small  H. 

The  proof  of  the  other  lower  bound  is  quite  similar 

Proof  of  Theorem  1:  The  bounds  on  the  norms  of  the  operators  follow 
immediately  from  Lemmzts  4  and  5. 

To  obtain  the  lower  bounds,  we  first  consider 

«=0 

Using  Lemma  5,  we  see  that  it  suffices  to  show  that 

t=0 

can  be  estimated  by 

CHA{u'*,u'*). 

By  the  definition  of  the  quadratic  forms 

A{u^  -  Q,u^,  Qiu^)     =     B{u^  -  Qiu'',  Q.u'") 
-S{u^-Q.u\Q,u^)    -    (c(ti''-Qiu''),Q.u''). 

By  using  the  definition  of  Qi,  the  first  term  of  the  right  hand  side  is  seen  to 
vanish. 

For  i  =  0,  the  second  term  can  be  estimated  by  CH'^''A{u'',u'*)  using  in- 
equality (iii)  and  Lemma  2.  We  note  that  S(Q,tz'',  Q.u'*)  =  0.  There  remains  to 
consider  5'(u'',  ^f^Q.u'').  By  using  the  inequality  (iii), 

I  X:  SCu"  -  Q,u\  Q.u")!  <  C\\u^\\a\\  E  ^.-^'IIi^- 

i=l  »=1 

Since,  for  each  a:  €  fi,  the  number  of  terms  Q.-tz''  that  differ  from  zero  is  imiformly 
bounded,  the  second  factor  can  be  estimated  by  {J2^i  l|Q«"''llis)^^^-  By  an  ele- 
mentary estimate,  which  shows  that  the  smallest  eigenvalue  of  Dirichlet  problem 
for  —A  on  Q,[  is  on  the  order  of  H^^,  and  Lemma  5,  the  required  inequality  is 
established. 

The  third  term  is  written  as  the  difference  of  two  expressions,  which  can  be 
handled  by  exactly  the  same  tools. 

The  estimate  for  the  operator  Q^'^^  is  obtained  similarly. 
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Figure  1:  An  extended  subregion 

5      Numerical  Results 

In  this  section,  we  present  some  numerical  results  to  demonstrate  the  behavior  of 
our  additive  Schwaxz  algorithms  for  both  symmetric  and  nonsymmetric  indefinite 
boundary  value  problems  in  B? .  Numerical  results  for  positive  definite  problems, 
both  symmetric  and  nonsymmetric,  have  previously  been  given  in  [3], [4], [12]. 
We  consider  the  problem 

r  Lu    =    /    in     n   =   [0,l]x  [0,1], 
1^     u    =    0    on    5f2. 

The  right  hand  side  /  is  always  chosen  so  that  the  exact  solution  is  u  = 
xt'^ sin{rx)sin{ry).  The  coefficients  of  L  are  specified  later  for  each  problem. 

We  use  a  two-level  subdivision  of  fi  as  described  in  section  2.  The  subregion 
fi^  is  obtained  by  enlarging  the  triangle  17,  as  in  Figure  1.  In  this  extension,  the 
same  number  ovlp  of  /i-level  triangles  are  added  in  all  directions. 

In  our  experiments  all  the  subproblems  are  solved  exactly  by  using  a  band 
solver  from  UNPACK.  We  stop  the  GMRES  method  as  soon  as  ||r,|U/||ro|U  < 
10~^.  We  have  foimd  that  the  overall  error  is  not  substantially  reduced  by  a 
more  stringent  stopping  criterion.  In  oui  tables,  the  error  denotes  the  difference 
between  the  computed  solution  and  the  exact  solution  of  the  continuous  problem 
measured  in  the  norms  indicated.  The  programs  have  been  run  in  single  precision 
on  the  Multiflow  computer  at  Yale  University. 
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Iteration 

A  —  norm;  residual 

L^  norm;  error 

L°^  norm;  error 

1 

2.50722 

0.318660 

0.719067 

2 

1.44028 

0.182950 

0.405074 

3 

0.971708 

1.63224E-02 

4.73967E-02 

4 

0.218693 

7.38034E-03 

2.46119E-02 

5 

5.54836E-02 

6.13811E-03 

1.94588E-02 

6 

3.40790E-02 

5.16731E-03 

1.53261E-02 

7 

2.60738E-02 

3.71766E-03 

9.73493E-03 

8 

1.87841E-02 

2.19535E-03 

6.23578E-03 

9 

1.03642E-02 

1.67804E-03 

4.90165E-03 

10 

6.81844E-03 

1.36607E-03 

4.00215E-03 

11 

5.02644E-03 

8.28570E-04 

2.39784E-03 

Table  1:   Convergence  history  for  Algorithm  1  and  Example  1.   Here  h   ^  =  75, 
H-'^  =  15,  ovlp  =  2  and  (5  =  IC.Ott^ 


Example  1.  We  consider  the  symmetric  and  indefinite  Helmholtz  equation 

(6) 


—  Au  —  6u    =    /    in     Cl 
u    =    0     on    do,, 


6  is  a  constant.  The  eigenvalues  of  the  operator  in  (6)  are  {P  +  P)t^^  —  S,  where 
i,j  are  positive  integers.  The  numerical  results  are  given  in  Tables  1  and  2. 


Example  2.  We  consider  a  nonsymmetric  and  indefinite  problem 
f  —  A  u  —  r}{du/dx  +  du/dy)  —  Su    =    f    in     U 

\  «  =  0  on  an, 


(7) 


The  numerical  results  are  given  in  Tables  3  and  4. 

We  note  that  in  some  of  the  experiments,  the  rate  of  convergence  is  unsat- 
isfactory, but  that  the  rate  of  convergence  improves  considerably  by  decreasing 


6      Two  Other  Methods 

We  conclude  by  outlining  how  some  other  results,  previously  analyzed  for  the 
positive  definite,  symmetric  case,  can  be  extended  to  the  class  of  elliptic  problems 
considered  in  this  paper.  We  confine  our  discussion  to  problems  in  the  plane;  both 
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Case  # 

6 

h-' 

H-' 

ovlp 

Algorithm  1 

Algorithm  2 

1 

Stt^ 

15 

3 

2 

11 

12 

2 

30 

3 

4 

11 

12 

3 

45 

3 

6 

12 

12 

4 

60 

3 

8 

12 

12 

5 

15 

5 

1 

10 

10 

6 

30 

5 

2 

12 

12 

( 

45 

5 

3 

12 

12 

8 

60 

5 

4 

12 

12 

10 

167r2 

45 

15 

1 

10 

10 

11 

60 

15 

1 

11 

11 

12 

75 

15 

2 

11 

11 

13 

60 

5 

4 

44 

33 

14 

60 

10 

2 

17 

17 

15 

60 

20 

1 

8 

8 

16 

30^2 

60 

20 

1 

16 

16 

17 

80 

20 

1 

17 

18 

Table  2:  Example  1.  The  last  two  columns  give  the  number  of  GMRES  iterations. 

of  the  algorithms  considered  here  need  to  be  modified  considerably  in  order  to 
obtain  fast  methods  for  problems  in  three  dimensions. 

We  first  consider  a  basic  iterative  substructuring  method  for  problems  in  two 
dimensions;  cf.  Dryja  and  Widlund  [8],  [9].  For  problems  that  are  nonsjinmetric, 
but  positive  definite,  the  result  to  be  formulated  has  pre\'iously  been  obtained  by 
Cai  [3],[4]._ 

When  iterative  substructuring  methods  are  tised,  the  region  is  divided  into 
substructures  and  elements  as  in  Section  2.  While  originally  derived  differently, 
it  has  been  demonstrated  by  Dryja  and  Widlund  [8]  that  these  methods  can  be 
viewed  as  additive  Schwarz  methods.  We  adopt  this  point  of  view. 

In  defining  the  partition  of  the  finite  element  space  into  subspaces,  we  use 
the  same  coarse  space  V^  introduced  in  Section  4.  We  also  use  subspaces  corre- 
sponding to  the  subregions  f2,j  =  n.UrijU^j-  These  subregions  play  the  same 
role  as  the  Vl[  in  Section  4.  Here  f],  and  Vlj  are  adjacent  substructures  with  the 
common  edge  ?,_,.  The  local  subspaces  are  thus  V.^  =  Hl{?lij)  fl  V^. 

Compared  with  the  case  considered  previously,  we  use  less  overlap  in  the 
sense  that  only  the  elements  of  V^  can  differ  from  zero  at  the  vertices  of  the 
substructures.  This  is  reflected  in  a  poorer  bound  for  the  constant  of  Lemma  1, 

Cl  <  c(mst.{l  +  log{Hlh)f; 
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Iteration 

A  —  norm;  residual 

L^  norm;  error 

L°°  norm;  error 

1 

2.99430 

0.309833 

0.696816 

2 

1.82397 

0.234651 

0.529782 

3 

1.34039 

0.136604 

0.313758 

4 

0.905845 

7.27307E-02 

0.172788 

5 

0.585598 

4.46239E-02 

0.108047 

6 

0.407054 

2.78872E-02 

6.71409E-02 

7 

0.288880 

1.37858E-02 

3.38832E-02 

8 

0.180577 

8.21095E-03 

2.06587E-02 

9 

0.129736 

4.44917E-03 

1.15359E-02 

10 

8.77408E-02 

2.19873E-03 

6.17071E-03 

11 

5.48599E-02 

1.18357E-03 

3.88315E-03 

12 

3.46359E-02 

7.18704E-04 

2.24515E-03 

13 

2.29454E-02 

4.35330E-04 

1.39198E-03 

14 

1.35519E-02 

2.90249E-04 

1.03428E-03 

15 

8.90639E-03 

2.18270E-04 

6.67672E-04 

16 

5.98530E-03 

1.90636E-04 

5.61312E-04 

17 

3.89341E-03 

1.72248E-04 

5.31457E-04 

Table  3:  Convergence  history  for  Algorithm  1  and  Example  2.  Here  h  ^  =  120, 
if-i  =  20,  ovlp  =  2,r)  =  16.07r  and  6  =  le.OTr^ 
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Case  # 

parameters 

/i-i 

H-' 

ovlp 

Algorithm  1 

Algorithm  2 

1 

6  =  37r2 

15 

5 

1 

13 

12 

2 

30 

5 

2 

17 

14 

3 

45 

5 

3 

18 

14 

4 

60 

5 

4 

18 

14 

5 

60 

6 

3 

16 

14 

6 

60 

10 

2 

12 

11 

7 

T]  =  167r 
6  =  167r2 

45 

15 

1 

17 

13 

8 

60 

15 

1 

18 

14 

9 

75 

15 

2 

25 

17 

10 

60 

20 

1 

13 

11 

11 

80 

20 

1 

14 

12 

12 

100 

20 

2 

18 

14 

13 

120 

20 

2 

17 

14 

14 

7/  =  30;r 
S  =  307r2 

60 

20 

1 

24 

16 

15 

120 

20 

2 

35 

19 

16 

75 

25 

1 

17 

13 

17 

100 

25 

1 

18 

14 

18 

120 

30 

1 

15 

13 

Table  4:  Example  2.  The  last  two  columns  give  the  number  of  GMRES  iterations. 
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cf.  Dryja  and  Widlxind  [8].  Lemma  1  is  modified  accordingly.  The  rest  of  the 
proof  carries  over  without  change.  In  Theorem  2,  we  use  the  notation  Q  = 
Qo  +  E  Q.:- 

Theorem  2  For  the  iterative  substructuring  method,  introduced  as  an  additive 
Schwarz  method  with  the  subspaces  V^  and  Vj^,  there  exist  constants  Hq  >  0, 
c{Ho)    >    0  and  C{Ho),  such  that  if  H  <  Hq, 

c{Ho){l  +  logiH/h))-''A{u\u'')     <     A{u^,Qu^) 

and 

A{Qu\Qu'')     <     C{Ho)A{u\u''). 

We  finally  show  that  the  result,  obtained  by  Yserentant  [24]  for  positive  defi- 
nite symmetric  problems,  can  be  extended  in  the  same  way.  We  note  that  Bank 
and  Yserentant  [1]  have  already  reported  on  successful  numerical  experiments 
with  an  accelerated  variant  of  this  algorithm  for  elliptic  problems  of  the  type 
considered  in  this  paper. 

We  assTime  that  the  region  fi  is  a  plane  polygon.  A  coarse  triangtdation  is 
introduced  as  before.  Its  triangles  are  recursively  subdivided  into  four  congruent 
triangles,  a  total  of  j  times.  The  characteristic  mesh  size  for  the  level  k  trian- 
gualtion  is  hk-  As  demonstrated  in  Yserentant  [8],  more  complicated  situations 
can  also  be  considered,  where  the  final  triangulation  is  highly  nonuniform,  but 
to  simplify  our  discussion,  we  only  consider  the  regular  case  in  this  paper. 

As  shown  in  Dryja  and  Widlund  [10],  Yserentant's  method  can  also  be  viewed 
as  an  additive  Schwarz  method  defined  by  a  set  of  subspaces.  Let  IkV  =  Ihi,v  be 
the  linear  interpolant  oi  v  E  V'*  onto  the  space  of  finite  elements  on  the  level  k 
triangulation.  The  following  identity  holds 

v  =  IoV-\-  {hv  -  lov)  +  ■■•  +  {IjV  -  Ij.iv)  ,  Vu  G  V''  . 

We  represent  V^  as 

where  Vq  =  V^  and,  for  A:  >  0,  V^  =  R{Ik  —  h-i)  is  the  range  of  the  operator 
(/t  —  Ik-i).  An  additive  Schwarz  method  is  defined  for  this  set  of  subspaces.  We 
obtain  Yserentant's  method  by,  for  A;  >  0,  replacing  the  resulting  problems  on 
the  subspaces  by  suitable  preconditioners. 

The  following  result  holds  for  the  family  of  elliptic  problems  considered  in 
this  paper.  Here  Q  denotes  the  operator  of  the  transformed  equation,  which 
corresponds  to  Yserentant's  method. 

Theorem  3  For  Yserentant's  method  there  exist  constants  Hq  >  0,  c{Ho)  >  0 
and  C{Ho),  such  that  if  H  <  Hq, 

c{Ho)j-^Aiu\u'')     <     A{u\Qu'') 
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and 

A{Qu\Qu'')     <     C(/fo)A(u^u'•). 

We  will  only  outline  how  this  result  can  be  established.  We  model  our  proof 
on  that  of  Theorem  1.  We  have  to  show  that  the  different  lemmas  hold  for  the 
spaces  just  introduced.  It  is  shown  in  Yserentant  [24]  that  Lemma  1  holds  with 
Co  <  P;  cf.  Lemmas  2.4  and  2.5  of  [24].  Lemma  2  is  still  valid  since  the  same 
coarse  operator  is  used  in  all  the  methods  considered  in  this  paper.  A  counter 
part  of  Lemma  3  can  be  obtained  as  well,  by  using  Lemma  2.4  of  [24].  Lemma  4 
is  modified  by  using  Lemma  2.7  of  [24],  a  result  that  makes  it  possible  to  obtain 
a  sharp  upper  bound  in  Yserentant's  main  theorem.  Lemma  5  can  be  modified 
straightforwardly.  One  change  is  required  in  the  proof  of  the  theorem.  The 
factor  II  E!Ii  Qiw''||L2  must  be  estimated  differently,  since,  typically,  all  these 
terms  differ  from  zero  everywhere.  By  using  Yserentant's  tools,  it  is  however 
possible  to  show  that 

WQ.AIl^    <   Ch,\\Q,u''\\A  . 

Since  the  h,  decay  geometrically,  the  triangle  and  Cauchy-Schwarz  inequalities 
give 


\\{:q.a\i  <  cH'j:\\Q.u'rA, 

i=l  1=1 

and  the  proof  can  be  completed. 
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